one non-mover

AFNI approach

r1 <- read.delim("~/Downloads/comp/3dmotion_run-01.1D", header = F, sep = "") %>% 
  mutate(run = "run-01")
r2 <- read.delim("~/Downloads/comp/3dmotion_run-02.1D", header = F, sep = "")%>% 
  mutate(run = "run-02")
# The output is in 9 ASCII formatted columns:
#   n  roll  pitch  yaw  dS  dL  dP  rmsold rmsnew

r_afni_p <- rbind(r1,r2) %>% 
  select(-V1) %>% 
  mutate(euc_dist = sqrt(V2^2+V3^2+V4^2+V5^2+V6^2+V7^2))
names(r_afni_p)[c(1:6,10)] <- paste0("afni_",c("x","y","z","roll","pitch","yaw","dist"))

r_afni<- rbind(r1,r2) %>% 
  group_by(run) %>% 
  mutate(euc_dist = sqrt(V2^2+V3^2+V4^2+V5^2+V6^2+V7^2),
         tr = row_number())

ggplot(r_afni, aes(x = tr, y = euc_dist, color = run))+
  geom_line()

r_afni <- r_afni%>% 
  group_by(run) %>% 
  mutate(euc_dist = sqrt(V2^2+V3^2+V4^2+V5^2+V6^2+V7^2),
         tr = row_number()) %>% 
  mutate_at(vars(V2:V7),~ .x - lag(.x))%>%
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>% 
  mutate(enorm = sqrt(V2_derivative^2+V3_derivative^2+V4_derivative^2+
                        V5_derivative^2+V6_derivative^2+V7_derivative^2))



print(paste("the mean Euclidean distance of run 1 using afni is: ",mean(r_afni$enorm[r_afni$run == "run-01"],na.rm = T)))
## [1] "the mean Euclidean distance of run 1 using afni is:  0.101661318128867"
print(paste("the mean Euclidean distance of run 2 using afni is: ",mean(r_afni$enorm[r_afni$run == "run-02"],na.rm = T)))
## [1] "the mean Euclidean distance of run 2 using afni is:  0.0850977083319949"

FreeSurfer approach

r1 <- read.delim("/Users/yanyan/Downloads/freesurfer_justin/functional/bn-run01/bold/003/fmcpr.mcdat", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-01")
r2 <- read.delim("/Users/yanyan/Downloads/freesurfer_justin/functional/bn-run02/bold/004/fmcpr.mcdat", header = F, sep = "")%>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-02")

r_fs <- rbind(r1, r2) %>% 
  group_by(run) %>% 
  # Re-reference V2-V7 to the first row in each run
  mutate(across(V2:V7, ~ .x - first(.x))) %>% 
  mutate(euc_dist = sqrt(V2^2 + V3^2 + V4^2 + V5^2 + V6^2 + V7^2),tr = row_number()) %>% 
  mutate(euc_dist = euc_dist - min(euc_dist)) %>% 
  mutate_at(vars(V2:V7), ~ .x - lag(.x)) %>% 
  mutate(across(V2:V7, ~ replace_na(.x, 0))) %>% 
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>% 
  mutate(enorm = sqrt(V2_derivative^2 + V3_derivative^2 + V4_derivative^2 +
      V5_derivative^2 + V6_derivative^2 + V7_derivative^2))

r_fs_p<- rbind(r1,r2) %>% 
  select(-V1) %>% 
  group_by(run) %>% 
  mutate(across(V2:V7, ~ .x - first(.x))) %>% 
  mutate(euc_dist = sqrt(V2^2 + V3^2 + V4^2 + V5^2 + V6^2 + V7^2),tr = row_number())

names(r_fs_p)[c(1:6,11)] <- paste0("fs_",c("x","y","z","roll","pitch","yaw","dist"))
print(paste("the mean Euclidean norm of run 1 using Freesurfer is: ",mean(r_fs$enorm[r_fs$run == "run-01"],na.rm = T)))
## [1] "the mean Euclidean norm of run 1 using Freesurfer is:  0.0783179186035793"
print(paste("the mean Euclidean norm of run 2 using Freesurfer is: ",mean(r_fs$enorm[r_fs$run == "run-02"],na.rm = T)))
## [1] "the mean Euclidean norm of run 2 using Freesurfer is:  0.0656547200134511"

#FSL-mutual info

r1 <- read.delim("/Users/yanyan/Desktop/MVPA/for_justin/pilot_bn230418/run-01/mid_run-01_mc_mutualinfo.nii.gz.par", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-01") %>% 
  mutate(tr = row_number())%>% 
  mutate_at(vars(V1:V3), ~ .x * 100)
r2 <- read.delim("/Users/yanyan/Desktop/MVPA/for_justin/pilot_bn230418/run-02/mid_run-02_mc_mutualinfo.nii.gz.par", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-02") %>% 
  mutate(tr = row_number())%>% 
  mutate_at(vars(V1:V3), ~ .x * 100)

r_fsl_mi_p<- rbind(r1,r2) %>% 
  group_by(run) %>% 
  mutate(across(V1:V6, ~ .x - first(.x))) %>% 
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2))
names(r_fsl_mi_p)[c(1:6,9)] <- paste0("fslmi_",c("y","x","z","roll","pitch","yaw","dist"))

r_fsl_mi <- rbind(r1,r2) %>%
  group_by(run) %>% 
  mutate(across(V1:V6, ~ .x - first(.x)),
         tr = row_number()) %>%  # Re-reference V2-V7 to the first row in each run
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2)) %>% 
  mutate(euc_dist = euc_dist - min(euc_dist)) %>% 
  mutate_at(vars(V1:V6), ~ .x - lag(.x)) %>% 
  mutate(across(V1:V6, ~ replace_na(.x, 0))) %>% 
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>%
  mutate(enorm = sqrt(V1_derivative^2+V2_derivative^2+V3_derivative^2+V4_derivative^2+
                        V5_derivative^2+V6_derivative^2))

#FSL-norm corr

r1 <- read.delim("/Users/yanyan/Desktop/MVPA/for_justin/pilot_bn230418/run-01/mid_run-01_mc_normcorr.nii.gz.par", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-01") %>% 
  mutate(tr = row_number()) %>% 
  mutate_at(vars(V1:V3), ~ .x * 100)
r2 <- read.delim("/Users/yanyan/Desktop/MVPA/for_justin/pilot_bn230418/run-02/mid_run-02_mc_normcorr.nii.gz.par", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-02") %>% 
  mutate(tr = row_number())%>% 
  mutate_at(vars(V1:V3), ~ .x * 100)

r_fsl_nc_p<- rbind(r1,r2) %>% 
  group_by(run) %>%
  mutate(across(V1:V6, ~ .x - first(.x))) %>% 
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2))
names(r_fsl_nc_p)[c(1:6,9)] <- paste0("fslnc_",c("y","x","z","roll","pitch","yaw","dist"))

r_fsl_nc <- rbind(r1,r2) %>%
  group_by(run) %>% 
  mutate(across(V1:V6, ~ .x - first(.x)),
         tr = row_number()) %>%  # Re-reference V2-V7 to the first row in each run
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2)) %>% 
  mutate(euc_dist = euc_dist - min(euc_dist)) %>% 
  mutate_at(vars(V1:V6), ~ .x - lag(.x)) %>% 
  mutate(across(V1:V6, ~ replace_na(.x, 0))) %>% 
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>%
  mutate(enorm = sqrt(V1_derivative^2+V2_derivative^2+V3_derivative^2+V4_derivative^2+
                        V5_derivative^2+V6_derivative^2))

fmriprep

r1 <- read.delim("~/Desktop/midmvpa/mid_data/derivatives/sub-bn230418/func/sub-bn230418_task-mid_run-1_desc-confounds_timeseries.tsv", header = T, sep = "\t") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-01")%>% 
  rowwise() %>% 
  mutate(motion_censor = as.integer(any(c_across(starts_with("motion_outlier")) == 1))) %>%
  select(global_signal:rmsd,trans_x:rot_z_power2,run:motion_censor)
## Warning: There were 22 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `global_signal_derivative1 = (structure(function (..., .x = ..1,
##   .y = ..2, . = ..1) ...`.
## Caused by warning in `as.numeric()`:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 21 remaining warnings.
r2 <- read.delim("~/Desktop/midmvpa/mid_data/derivatives/sub-bn230418/func/sub-bn230418_task-mid_run-2_desc-confounds_timeseries.tsv", header = T, sep = "\t")%>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-02")%>% 
  rowwise() %>% 
  mutate(motion_censor = as.integer(any(c_across(starts_with("motion_outlier")) == 1))) %>%
  select(global_signal:rmsd,trans_x:rot_z_power2,run:motion_censor)
## Warning: There were 22 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `global_signal_derivative1 = (structure(function (..., .x = ..1,
##   .y = ..2, . = ..1) ...`.
## Caused by warning in `as.numeric()`:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 21 remaining warnings.
r_fp_p<- rbind(r1,r2) %>% 
  select(trans_x,trans_y,trans_z, rot_x,rot_y,rot_z) %>% 
  mutate_at(vars(rot_x,rot_y,rot_z),~ 50 * .x) %>% 
  mutate(euc_dist = sqrt(trans_x^2+trans_y^2+trans_z^2+rot_x^2+rot_y^2+rot_z^2))

names(r_fp_p)[1:7] <- paste0("fp_",c("roll","pitch","yaw","x","y","z","dist"))

r_fp<- rbind(r1,r2) %>% 
  group_by(run) %>% 
  mutate(#enorm = framewise_displacement,  
         tr = row_number()) %>% 
  mutate_at(vars(rot_x,rot_y,rot_z),~ 50 * .x) %>% 
  rename(V1=trans_x,V2=trans_y,V3=trans_z, V4=rot_x,V5=rot_y,V6=rot_z) %>% 
  group_by(run) %>% 
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2),
         tr = row_number()) %>% 
  mutate(euc_dist = euc_dist-min(euc_dist))%>% 
  mutate_at(vars(V1:V6),~ .x - lag(.x))%>%
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>% 
  mutate(enorm = sqrt(V1_derivative^2+V2_derivative^2+V3_derivative^2+V4_derivative^2+
                        V5_derivative^2+V6_derivative^2))

ggplot(r_fp, aes(x = framewise_displacement, y = enorm, color = run))+
  geom_point()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(r_fp, aes(x = tr, y = euc_dist, color = run))+
  geom_line()

ggplot(r_fp, aes(x = tr, y = enorm, color = run))+
  geom_line()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

print(paste("the mean Euclidean norm of run 1 using Freesurfer is: ",mean(r_fp$enorm[r_fp$run == "run-01"],na.rm = T)))
## [1] "the mean Euclidean norm of run 1 using Freesurfer is:  0.116737959117432"
print(paste("the mean Euclidean norm of run 2 using Freesurfer is: ",mean(r_fp$enorm[r_fp$run == "run-02"],na.rm = T)))
## [1] "the mean Euclidean norm of run 2 using Freesurfer is:  0.0742764342741759"

plot all time series

ggplot(r_afni, aes(x = tr, y = euc_dist, color = run))+
  geom_line()+
  ylim(0,1)+
  labs(title = "AFNI")+
  
  ggplot(r_fsl_mi, aes(x = tr, y = euc_dist, color = run))+
  labs(title = "FSL coreg - mutualinfo")+
  geom_line()+
  
  ggplot(r_fsl_nc, aes(x = tr, y = euc_dist, color = run))+
  labs(title = "FSL coreg - normcorr")+
  geom_line()+
  
  ggplot(r_fs, aes(x = tr, y = euc_dist, color = run))+
  labs(title = "FreeSurfer coreg")+
  geom_line()+
  
    ggplot(r_fp, aes(x = tr, y = euc_dist, color = run))+
  labs(title = "fmriPrep coreg")+
  geom_line()+
  ggplot(r_afni, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+
  

  ggplot(r_fsl_mi, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+

  ggplot(r_fsl_nc, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+
  
  ggplot(r_fs, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+

  ggplot(r_fp, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+
  plot_layout(nrow = 2)
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

compare estimates

params <- cbind(r_fs_p,
           r_afni_p[,c(1:6,10)],
           r_fsl_mi_p[,c(1:6,9)],
           r_fsl_nc_p[,c(1:6,9)],
           r_fp_p[,1:7]) %>% 
  ungroup() %>% 
  mutate(tr = row_number()) %>% 
  select(-V8:-run) %>% 
  relocate(tr) %>% 
  # mutate_at(vars(fs_x:fp_yaw), ~ . - lag(.)) %>% 
  pivot_longer(fs_x:fp_dist,names_to = c("method","param"), values_to = "value", names_sep = "_")

ggplot(params %>% filter(param %in% c("dist")))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y")

ggplot(params %>% filter(param %in% c("x","y","z")))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y")

ggplot(params %>% filter(param %in% c("roll","pitch","yaw")))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y")+
  scale_color_brewer(palette = "Set2")

params_ <- cbind(r_fs_p,
           r_afni_p[,c(1:6,10)],
           r_fsl_mi_p[,c(1:6,9)],
           r_fsl_nc_p[,c(1:6,9)],
           r_fp_p[,1:7]) %>% 
  ungroup() %>% 
  mutate(tr = row_number()) %>% 
  select(-V8:-run) %>% 
  relocate(tr) %>% 
  mutate_at(vars(fs_x:fp_dist), ~ . - lag(.)) %>% 
  mutate(across(fs_x:fp_dist, ~ replace_na(.x, 0))) %>% 
  pivot_longer(fs_x:fp_dist,names_to = c("method","param"), values_to = "value", names_sep = "_")

ggplot(params_ %>% filter(param %in% c("dist"), tr != 2, tr != 258))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y", ncol = 1)

ggplot(params_ %>% filter(param %in% c("x","y","z"), tr != 2, tr != 258))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y", ncol = 1)

ggplot(params_ %>% filter(param %in% c("roll","pitch","yaw"), tr != 2, tr != 258))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y", ncol = 1)+
  scale_color_brewer(palette = "Set2")

d1<-params_ %>% 
  filter(param == "dist") %>% 
         pivot_wider(names_from = c("method","param"), values_from = c("value")) %>% 
  filter(tr != 1) %>% 
  select(-tr)

PerformanceAnalytics::chart.Correlation(d1)
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter

corr <- round(cor(d1, use = "complete.obs"), 1)
p.mat <- rstatix::cor_pmat(d1)

ggcorrplot(
  corr, 
  hc.order = TRUE,
  type = "lower",
  outline.color = "white",
  ggtheme = ggplot2::theme_gray,lab_size = 4,
  colors = c("#6D9EC1", "white", "#E46726"),
  lab = TRUE,
  title = "Euclidean Distance"
)

plot all enorms

r <- cbind(r_fs %>%  ungroup() %>% select(enorm) %>% rename(enorm_fs = enorm),
           r_afni %>%  ungroup()%>% select(enorm) %>% rename(enorm_afni = enorm),
           r_fsl_mi %>%  ungroup()%>% select(enorm) %>% rename(enorm_fsl_mi = enorm),
           r_fsl_nc %>%  ungroup()%>% select(enorm) %>% rename(enorm_fsl_nc = enorm),
           r_fp %>%  ungroup()%>% select(enorm) %>% rename(enorm_fp = enorm)) %>% 
  mutate(tr = row_number())


(ggplot(r)+
    geom_line(aes(x = tr, y = enorm_afni), color = 'orange')+
    geom_line(aes(x = tr, y = enorm_fs), color = 'skyblue')+
    geom_line(aes(x = tr, y = enorm_fsl_mi), color = 'blue4')+
    geom_line(aes(x = tr, y = enorm_fsl_nc), color = 'purple')+
    geom_line(aes(x = tr, y = enorm_fp), color = 'salmon')+
    annotate("text",x = 0, y = 0.3, color = 'orange', label = "AFNI")+
    annotate("text",x = 0, y = 0.33, color = 'skyblue', label = "FreeSurfer")+
    annotate("text",x = 0, y = 0.27, color = 'blue4', label = "FSL-mutual")+
    annotate("text",x = 0, y = 0.24, color = 'purple', label = "FSL-normcorr")+
    annotate("text",x = 0, y = 0.21, color = 'salmon', label = "fmriPrep")+
    labs(y = "eucl. norm",
         title = "bn230418"))/
  
  (  ggplot(r, aes(x = enorm_afni, y = enorm_fs))+
       geom_point(size = 0.1)+
       geom_smooth(method = "lm")|
       
       ggplot(r, aes(x = enorm_afni, y = enorm_fsl_mi))+
       geom_point(size = 0.1)+
       geom_smooth(method = "lm")|
       
       ggplot(r, aes(x = enorm_fs, y = enorm_fsl_mi))+
       geom_point(size = 0.1)+
       geom_smooth(method = "lm")|
       
       ggplot(r, aes(x = enorm_afni, y = enorm_fp))+
       geom_point(size = 0.1)+
       geom_smooth(method = "lm"))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

corr <- round(cor(r %>% select(enorm_fs:enorm_fp), use = "complete.obs"), 1)
p.mat <- rstatix::cor_pmat(r %>% select(enorm_fs:enorm_fp))

ggcorrplot(
  corr, 
  hc.order = TRUE,
  type = "lower",
  outline.color = "white",
  ggtheme = ggplot2::theme_gray,lab_size = 4,
  colors = c("#6D9EC1", "white", "#E46726"),
  lab = TRUE
)

one mover

AFNI approach

r1 <- read.delim("~/Desktop/MVPA/for_justin/cw231123/3dmotion_run-01.1D", header = F, sep = "") %>% 
  mutate(run = "run-01")
r2 <- read.delim("~/Desktop/MVPA/for_justin/cw231123/3dmotion_run-02.1D", header = F, sep = "")%>% 
  mutate(run = "run-02")
# The output is in 9 ASCII formatted columns:
#   n  roll  pitch  yaw  dS  dL  dP  rmsold rmsnew

r_afni_p <- rbind(r1,r2) %>% 
  select(-V1) %>% 
  mutate(euc_dist = sqrt(V2^2+V3^2+V4^2+V5^2+V6^2+V7^2))
names(r_afni_p)[c(1:6,10)] <- paste0("afni_",c("x","y","z","roll","pitch","yaw","dist"))

r_afni<- rbind(r1,r2) %>% 
  group_by(run) %>% 
  mutate(euc_dist = sqrt(V2^2+V3^2+V4^2+V5^2+V6^2+V7^2),
         tr = row_number())

ggplot(r_afni, aes(x = tr, y = euc_dist, color = run))+
  geom_line()

r_afni <- r_afni%>% 
  group_by(run) %>% 
  mutate(euc_dist = sqrt(V2^2+V3^2+V4^2+V5^2+V6^2+V7^2),
         tr = row_number()) %>% 
  mutate_at(vars(V2:V7),~ .x - lag(.x))%>%
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>% 
  mutate(enorm = sqrt(V2_derivative^2+V3_derivative^2+V4_derivative^2+
                        V5_derivative^2+V6_derivative^2+V7_derivative^2))

FreeSurfer approach

r1 <- read.delim("/Users/yanyan/Downloads/freesurfer_justin/functional/cw-run01/bold/001/fmcpr.mcdat", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-01")
r2 <- read.delim("/Users/yanyan/Downloads/freesurfer_justin/functional/cw-run02/bold/002/fmcpr.mcdat", header = F, sep = "")%>% 
  mutate_all(~as.numeric(.x)) %>% 
   mutate(run = "run-02")


r_fs <- rbind(r1, r2) %>% 
  group_by(run) %>% 
  # Re-reference V2-V7 to the first row in each run
  mutate(across(V2:V7, ~ .x - first(.x))) %>% 
  mutate(euc_dist = sqrt(V2^2 + V3^2 + V4^2 + V5^2 + V6^2 + V7^2),tr = row_number()) %>% 
  mutate(euc_dist = euc_dist - min(euc_dist)) %>% 
  mutate_at(vars(V2:V7), ~ .x - lag(.x)) %>% 
  mutate(across(V2:V7, ~ replace_na(.x, 0))) %>% 
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>% 
  mutate(enorm = sqrt(V2_derivative^2 + V3_derivative^2 + V4_derivative^2 +
      V5_derivative^2 + V6_derivative^2 + V7_derivative^2))

r_fs_p<- rbind(r1,r2) %>% 
  select(-V1) %>% 
  group_by(run) %>% 
  mutate(euc_dist = sqrt(V2^2 + V3^2 + V4^2 + V5^2 + V6^2 + V7^2),tr = row_number()) %>% 
  mutate(euc_dist = euc_dist - min(euc_dist))

names(r_fs_p)[c(1:6,11)] <- paste0("fs_",c("x","y","z","roll","pitch","yaw","dist"))

FSL-mutual info

r1 <- read.delim("/Users/yanyan/Desktop/MVPA/for_justin/cw231123/run-01/mid_run-01_mc_mutualinfo.nii.gz.par", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-01") %>% 
  mutate(rowid = row_number())%>% 
  mutate_at(vars(V1:V3), ~ .x * 100)
r2 <- read.delim("/Users/yanyan/Desktop/MVPA/for_justin/cw231123/run-02/mid_run-02_mc_mutualinfo.nii.gz.par", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-02") %>% 
  mutate(rowid = row_number())%>% 
  mutate_at(vars(V1:V3), ~ .x * 100)

r_fsl_mi_p<- rbind(r1,r2) %>% 
  group_by(run) %>% 
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2))
names(r_fsl_mi_p)[c(1:6,9)] <- paste0("fslmi_",c("y","x","z","roll","pitch","yaw","dist"))

r_fsl_mi <- rbind(r1,r2) %>%
  group_by(run) %>% 
  mutate(across(V1:V6, ~ .x - first(.x)),
         tr = row_number()) %>%  # Re-reference V2-V7 to the first row in each run
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2)) %>% 
  mutate(euc_dist = euc_dist - min(euc_dist)) %>% 
  mutate_at(vars(V1:V6), ~ .x - lag(.x)) %>% 
  mutate(across(V1:V6, ~ replace_na(.x, 0))) %>% 
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>%
  mutate(enorm = sqrt(V1_derivative^2+V2_derivative^2+V3_derivative^2+V4_derivative^2+
                        V5_derivative^2+V6_derivative^2))

FSL norm corr

r1 <- read.delim("/Users/yanyan/Desktop/MVPA/for_justin/cw231123/run-01/mid_run-01_mc_normcorr.nii.gz.par", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-01") %>% 
  mutate(rowid = row_number())%>% 
  mutate_at(vars(V1:V3), ~ .x * 100)
r2 <- read.delim("/Users/yanyan/Desktop/MVPA/for_justin/cw231123/run-02/mid_run-02_mc_normcorr.nii.gz.par", header = F, sep = "") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-02") %>% 
  mutate(rowid = row_number())%>% 
  mutate_at(vars(V1:V3), ~ .x * 100)


r_fsl_nc_p<- rbind(r1,r2) %>% 
  group_by(run) %>%
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2))
names(r_fsl_nc_p)[c(1:6,9)] <- paste0("fslnc_",c("y","x","z","roll","pitch","yaw","dist"))

r_fsl_nc <- rbind(r1,r2) %>%
  group_by(run) %>% 
  mutate(across(V1:V6, ~ .x - first(.x)),
         tr = row_number()) %>%  # Re-reference V2-V7 to the first row in each run
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2)) %>% 
  mutate(euc_dist = euc_dist - min(euc_dist)) %>% 
  mutate_at(vars(V1:V6), ~ .x - lag(.x)) %>% 
  mutate(across(V1:V6, ~ replace_na(.x, 0))) %>% 
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>%
  mutate(enorm = sqrt(V1_derivative^2+V2_derivative^2+V3_derivative^2+V4_derivative^2+
                        V5_derivative^2+V6_derivative^2))

fmriprep

r1 <- read.delim("~/Desktop/midmvpa/mid_data/derivatives/sub-cw231123/func/sub-cw231123_task-mid_run-1_desc-confounds_timeseries.tsv", header = T, sep = "\t") %>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-01")%>% 
  rowwise() %>% 
  mutate(motion_censor = as.integer(any(c_across(starts_with("motion_outlier")) == 1))) %>%
  select(global_signal:rmsd,trans_x:rot_z_power2,run:motion_censor)
## Warning: There were 22 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `global_signal_derivative1 = (structure(function (..., .x = ..1,
##   .y = ..2, . = ..1) ...`.
## Caused by warning in `as.numeric()`:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 21 remaining warnings.
r2 <- read.delim("~/Desktop/midmvpa/mid_data/derivatives/sub-cw231123/func/sub-cw231123_task-mid_run-2_desc-confounds_timeseries.tsv", header = T, sep = "\t")%>% 
  mutate_all(~as.numeric(.x)) %>% 
  mutate(run = "run-02")%>% 
  rowwise() %>% 
  mutate(motion_censor = as.integer(any(c_across(starts_with("motion_outlier")) == 1))) %>%
  select(global_signal:rmsd,trans_x:rot_z_power2,run:motion_censor)
## Warning: There were 22 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `global_signal_derivative1 = (structure(function (..., .x = ..1,
##   .y = ..2, . = ..1) ...`.
## Caused by warning in `as.numeric()`:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 21 remaining warnings.
r_fp_p<- rbind(r1,r2) %>% 
  select(trans_x,trans_y,trans_z, rot_x,rot_y,rot_z) %>% 
   mutate_at(vars(rot_x,rot_y,rot_z),~ 50 * .x) %>% 
  mutate(euc_dist = sqrt(trans_x^2+trans_y^2+trans_z^2+rot_x^2+rot_y^2+rot_z^2))

names(r_fp_p)[1:7] <- paste0("fp_",c("roll","pitch","yaw","x","y","z","dist"))

r_fp<- rbind(r1,r2) %>% 
  group_by(run) %>% 
  mutate(#enorm = framewise_displacement,  
    tr = row_number()) %>% 
   mutate_at(vars(rot_x,rot_y,rot_z),~ 50 * .x) %>% 
  rename(V1=trans_x,V2=trans_y,V3=trans_z, V4=rot_x,V5=rot_y,V6=rot_z) %>% 
  group_by(run) %>% 
  mutate(euc_dist = sqrt(V1^2+V2^2+V3^2+V4^2+V5^2+V6^2),
         tr = row_number()) %>% 
  mutate(euc_dist = euc_dist-min(euc_dist))%>% 
  mutate_at(vars(V1:V6),~ .x - lag(.x))%>%
  rename_with(~ paste0(.x, "_derivative"), matches("V")) %>% 
  mutate(enorm = sqrt(V1_derivative^2+V2_derivative^2+V3_derivative^2+V4_derivative^2+
                        V5_derivative^2+V6_derivative^2))

ggplot(r_fp, aes(x = framewise_displacement, y = enorm, color = run))+
  geom_point()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

ggplot(r_fp, aes(x = tr, y = euc_dist, color = run))+
  geom_line()

ggplot(r_fp, aes(x = tr, y = enorm, color = run))+
  geom_line()
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_line()`).

print(paste("the mean Euclidean norm of run 1 using Freesurfer is: ",mean(r_fp$enorm[r_fp$run == "run-01"],na.rm = T)))
## [1] "the mean Euclidean norm of run 1 using Freesurfer is:  0.1428224399266"
print(paste("the mean Euclidean norm of run 2 using Freesurfer is: ",mean(r_fp$enorm[r_fp$run == "run-02"],na.rm = T)))
## [1] "the mean Euclidean norm of run 2 using Freesurfer is:  0.100254907705082"

plot all time series

ggplot(r_afni, aes(x = tr, y = euc_dist, color = run))+
  geom_line()+
  ylim(0,1)+
  labs(title = "AFNI")+
  
  ggplot(r_fsl_mi, aes(x = tr, y = euc_dist, color = run))+
  labs(title = "FSL coreg - mutualinfo")+
  geom_line()+
  
  ggplot(r_fsl_nc, aes(x = tr, y = euc_dist, color = run))+
  labs(title = "FSL coreg - normcorr")+
  geom_line()+
  
  ggplot(r_fs, aes(x = tr, y = euc_dist, color = run))+
  labs(title = "FreeSurfer coreg")+
  geom_line()+
  
    ggplot(r_fp, aes(x = tr, y = euc_dist, color = run))+
  labs(title = "fmriPrep coreg")+
  geom_line()+
  ggplot(r_afni, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+
  

  ggplot(r_fsl_mi, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+

  ggplot(r_fsl_nc, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+
  
  ggplot(r_fs, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+

  ggplot(r_fp, aes(x = tr, y = enorm, color = run))+
  geom_line()+
  ylim(0,0.4)+
  plot_layout(nrow = 2)

compare estimates

params <- cbind(r_fs_p,
           r_afni_p[,c(1:6,10)],
           r_fsl_mi_p[,c(1:6,9)],
           r_fsl_nc_p[,c(1:6,9)],
           r_fp_p[,1:7]) %>% 
  ungroup() %>% 
  mutate(tr = row_number()) %>% 
  select(-V8:-run) %>% 
  relocate(tr) %>% 
  # mutate_at(vars(fs_x:fp_yaw), ~ . - lag(.)) %>% 
  pivot_longer(fs_x:fp_dist,names_to = c("method","param"), values_to = "value", names_sep = "_")

ggplot(params %>% filter(param %in% c("dist")))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y")

ggplot(params %>% filter(param %in% c("x","y","z")))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y")

ggplot(params %>% filter(param %in% c("roll","pitch","yaw")))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y")+
  scale_color_brewer(palette = "Set2")

params_ <- cbind(r_fs_p,
           r_afni_p[,c(1:6,10)],
           r_fsl_mi_p[,c(1:6,9)],
           r_fsl_nc_p[,c(1:6,9)],
           r_fp_p[,1:7]) %>% 
  ungroup() %>% 
  mutate(tr = row_number()) %>% 
  select(-V8:-run) %>% 
  relocate(tr) %>% 
  mutate_at(vars(fs_x:fp_dist), ~ . - lag(.)) %>% 
  mutate(across(fs_x:fp_dist, ~ replace_na(.x, 0))) %>% 
  pivot_longer(fs_x:fp_dist,names_to = c("method","param"), values_to = "value", names_sep = "_")

ggplot(params_ %>% filter(param %in% c("dist"), tr != 2, tr != 258))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y", ncol = 1)

ggplot(params_ %>% filter(param %in% c("x","y","z"), tr != 2, tr != 258))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y", ncol = 1)

ggplot(params_ %>% filter(param %in% c("roll","pitch","yaw"), tr != 2, tr != 258))+
  geom_line(aes(x = tr, y = value, color = param))+
  facet_wrap(~method, scales = "free_y", ncol = 1)+
  scale_color_brewer(palette = "Set2")

d1<-params_ %>% 
  filter(param == "dist") %>% 
         pivot_wider(names_from = c("method","param"), values_from = c("value")) %>% 
  filter(tr != 1) %>% 
  select(-tr)

PerformanceAnalytics::chart.Correlation(d1)
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter

corr <- round(cor(d1, use = "complete.obs"), 1)
p.mat <- rstatix::cor_pmat(d1)

ggcorrplot(
  corr, 
  hc.order = TRUE,
  type = "lower",
  outline.color = "white",
  ggtheme = ggplot2::theme_gray,lab_size = 4,
  colors = c("#6D9EC1", "white", "#E46726"),
  lab = TRUE,
  title = "Euclidean Distance"
)

d2<-params_ %>% 
  filter(param == "x") %>% 
         pivot_wider(names_from = c("method","param"), values_from = c("value")) %>% 
  filter(tr != 1) %>% 
  select(-tr)


corr <- round(cor(d2, use = "complete.obs"), 1)
p.mat <- rstatix::cor_pmat(d2)

ggcorrplot(
  corr, 
  hc.order = TRUE,
  type = "lower",
  outline.color = "white",
  ggtheme = ggplot2::theme_gray,lab_size = 4,
  colors = c("#6D9EC1", "white", "#E46726"),
  lab = TRUE,
  title = "X displacement"
)

plot all enorms

r <- cbind(r_fs %>%  ungroup() %>% select(enorm) %>% rename(enorm_fs = enorm),
           r_afni %>%  ungroup()%>% select(enorm) %>% rename(enorm_afni = enorm),
           r_fsl_mi %>%  ungroup()%>% select(enorm) %>% rename(enorm_fsl_mi = enorm),
           r_fsl_nc %>%  ungroup()%>% select(enorm) %>% rename(enorm_fsl_nc = enorm),
           r_fp %>%  ungroup()%>% select(enorm) %>% rename(enorm_fp = enorm)) %>% 
  mutate(tr = row_number())

(ggplot(r)+
    geom_line(aes(x = tr, y = enorm_afni), color = 'orange')+
    geom_line(aes(x = tr, y = enorm_fs), color = 'skyblue')+
    geom_line(aes(x = tr, y = enorm_fsl_mi), color = 'blue4')+
    geom_line(aes(x = tr, y = enorm_fsl_nc), color = 'purple')+
    geom_line(aes(x = tr, y = enorm_fp), color = 'salmon')+
    annotate("text",x = 0, y = 0.3, color = 'orange', label = "AFNI")+
    annotate("text",x = 0, y = 0.33, color = 'skyblue', label = "FreeSurfer")+
    annotate("text",x = 0, y = 0.27, color = 'blue4', label = "FSL-mutual")+
    annotate("text",x = 0, y = 0.24, color = 'purple', label = "FSL-normcorr")+
    annotate("text",x = 0, y = 0.21, color = 'salmon', label = "fmriPrep")+
    labs(y = "eucl. norm",
         title = "cw231123"))/
  
  (  ggplot(r, aes(x = enorm_afni, y = enorm_fs))+
       geom_point(size = 0.1)+
       geom_smooth(method = "lm")|
       
       ggplot(r, aes(x = enorm_afni, y = enorm_fsl_mi))+
       geom_point(size = 0.1)+
       geom_smooth(method = "lm")|
       
       ggplot(r, aes(x = enorm_fs, y = enorm_fsl_mi))+
       geom_point(size = 0.1)+
       geom_smooth(method = "lm")|
       
       ggplot(r, aes(x = enorm_afni, y = enorm_fp))+
       geom_point(size = 0.1)+
       geom_smooth(method = "lm"))
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## Removed 1 row containing missing values or values outside the scale range
## (`geom_line()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_smooth()`).
## Warning: Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).
## `geom_smooth()` using formula = 'y ~ x'
## `geom_smooth()` using formula = 'y ~ x'
## Warning: Removed 2 rows containing non-finite outside the scale range (`stat_smooth()`).
## Removed 2 rows containing missing values or values outside the scale range
## (`geom_point()`).

corr <- round(cor(r %>% select(enorm_fs:enorm_fp), use = "complete.obs"), 1)
p.mat <- rstatix::cor_pmat(r %>% select(enorm_fs:enorm_fp))

ggcorrplot(
  corr, 
  hc.order = TRUE,
  type = "lower",
  outline.color = "white",
  ggtheme = ggplot2::theme_gray,lab_size = 4,
  colors = c("#6D9EC1", "white", "#E46726"),
  lab = TRUE
)

compare the time course of nacc

read in data

my.methods = c("nomc","afni","fsl-mi","fsl-nc","freesurfer","fmriprep")
my.subjects = c("cw231123","bn230418")
r_all <- NULL
for (s in my.subjects){
  for (m in my.methods){
  r1 <- read.delim(paste0("~/Downloads/comp/mc_all/",m,"/roi_ts/",s,"_mid_b4_nacc8mm.1D"), header = F, sep = "") %>% 
    mutate(method = m,
           subject = s,
           tr = row_number())
  r_all <- rbind(r_all,r1)
  }
}

r_all <- r_all %>% 
  pivot_wider(names_from = "method", values_from = V1)

PerformanceAnalytics::chart.Correlation(r_all %>% filter(subject =="bn230418") %>% select(-subject,-tr))
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter

PerformanceAnalytics::chart.Correlation(r_all %>% filter(subject =="cw231123") %>% select(-subject,-tr))
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter
## Warning in par(usr): argument 1 does not name a graphical parameter

ggplot(r_all)+
  geom_line(aes(x = tr, y = nomc), color = 'orange', linewidth = 0.5)+
  geom_line(aes(x = tr, y = afni), color = 'skyblue', linewidth = 0.5)+
  geom_line(aes(x = tr, y = `fsl-mi`), color = 'purple', linewidth = 0.5)+
  geom_line(aes(x = tr, y = `fsl-nc`), color = 'violet', linewidth = 0.5)+
  geom_line(aes(x = tr, y = fmriprep), color = 'black', linewidth = 0.5)+

  annotate("text",x = -1, y = 2.5, color = 'orange', label = "no mc")+
  annotate("text",x = -1, y = 2, color = 'skyblue', label = "AFNI")+
  annotate("text",x = -1, y = 1, color = 'purple', label = "FSL-mutual")+
  annotate("text",x = -1, y = 1.5, color = 'black', label = "fmriprep")+
  annotate("text",x = -1, y = 0.5, color = 'violet', label = "FSL-normcorr")+
  theme_bw()+
  labs(y = "nacc % signal change",
       title = "effect of motion correction on nacc average signal",
       subtitle = "holding preprocessing constant")+
  facet_wrap(~subject,nrow = 2)

ggplot(r_all)+
  geom_line(aes(x = tr, y = nomc), color = 'orange', linewidth = 1.5)+
  geom_line(aes(x = tr, y = afni), color = 'skyblue', linewidth = 0.5)+
  annotate("text",x = -1, y = 2.5, color = 'orange', label = "no mc")+
  annotate("text",x = -1, y = 2, color = 'skyblue', label = "AFNI")+
  theme_bw()+
  facet_wrap(~subject,nrow = 2)

ggplot(r_all)+
  geom_line(aes(x = tr, y = nomc), color = 'orange', linewidth = 1.5)+
  geom_line(aes(x = tr, y = afni), color = 'skyblue', linewidth = 1)+
  geom_line(aes(x = tr, y = fmriprep), color = 'blue4', linewidth = 0.2)+
  annotate("text",x = -1, y = 2.5, color = 'orange', label = "no mc")+
  annotate("text",x = -1, y = 2, color = 'skyblue', label = "AFNI")+
  annotate("text",x = -1, y = 1.5, color = 'blue4', label = "fmriprep")+
  facet_wrap(~subject,nrow = 2, scales = "free_y")

ggplot(r_all)+
  geom_point(aes(x = afni, y = `fsl-mi`), color = 'orange')+
  facet_wrap(~subject,nrow = 2)

ggplot(r_all)+
  geom_point(aes(x = afni, y = fmriprep), color = 'orange')+
  facet_wrap(~subject,nrow = 2, scales = "free_y")

summary(lm(scale(fmriprep) ~ scale(afni),r_all %>% filter(subject == "bn230418")))
## 
## Call:
## lm(formula = scale(fmriprep) ~ scale(afni), data = r_all %>% 
##     filter(subject == "bn230418"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.6963 -0.6192  0.0053  0.6470  3.0339 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 8.763e-17  4.105e-02   0.000        1    
## scale(afni) 2.739e-01  4.108e-02   6.667 6.38e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9626 on 548 degrees of freedom
## Multiple R-squared:  0.07503,    Adjusted R-squared:  0.07334 
## F-statistic: 44.45 on 1 and 548 DF,  p-value: 6.381e-11
summary(lm(scale(fmriprep) ~ scale(afni),r_all %>% filter(subject == "cw231123")))
## 
## Call:
## lm(formula = scale(fmriprep) ~ scale(afni), data = r_all %>% 
##     filter(subject == "cw231123"))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6554 -0.5729 -0.0233  0.5873  2.8291 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.045e-18  3.926e-02   0.000        1    
## scale(afni)  3.921e-01  3.930e-02   9.979   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.9207 on 548 degrees of freedom
## Multiple R-squared:  0.1538, Adjusted R-squared:  0.1522 
## F-statistic: 99.59 on 1 and 548 DF,  p-value: < 2.2e-16
nacc_tc_b4 <- read_csv("~/Downloads/comp/nacc_tc/timecourses_b4_long.csv") %>% 
    filter(voi %in% c("nacc8mm")) %>%
  mutate(cue_value = factor(trialtype, levels = 1:6, labels = c("-$0","-$1","-$5","+$0","+$1","+$5"))) 
## Rows: 32400 Columns: 19
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## chr  (6): subject, method, cue_value, trial_gain, total, voi
## dbl (13): trial, trialonset, trialtype, target_ms, rt, hit, iti, drift, tota...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
nacc_tc_b4_all <- nacc_tc_b4 %>% 
  group_by(subject,voi,method) %>% 
  mutate(id = row_number())

ggplot(nacc_tc_b4_all %>% filter(voi == "nacc8mm"), aes(x = id, y = BOLD, group = method, 
                               fill = method))+
  geom_line(aes(color = method))+
  facet_wrap(~subject, nrow = 2)+
  theme_bw()+
  scale_fill_brewer(palette = "Paired")+
  scale_color_brewer(palette = "Paired")

nacc_tc_b4_summary <- nacc_tc_b4 %>% 
  filter(cue_value %in% c("+$0","+$5")) %>% 
  group_by(subject,voi,method,tr,cue_value) %>% 
  summarise(mBOLD = mean(BOLD,na.rm = T),
            seBOLD = sd(BOLD,na.rm = T)/ sqrt(n()))
## `summarise()` has grouped output by 'subject', 'voi', 'method', 'tr'. You can
## override using the `.groups` argument.
ggplot(nacc_tc_b4_summary, aes(x = tr, y = mBOLD, group = method, 
                               fill = method))+
  geom_line(aes(color = method,linetype = cue_value), linewidth = 1)+
  geom_ribbon(aes(group=method,ymin = mBOLD - seBOLD, ymax = mBOLD + seBOLD),
              alpha = 0.1)+
  facet_wrap(~subject+cue_value, nrow = 1)+
  theme_bw()+
  scale_x_continuous(breaks = 1:10)+
  labs(title = "average time course")

nacc_tc_b4_summary_cueval <- nacc_tc_b4 %>%
  group_by(subject,voi,method,tr,cue_value) %>%
  summarise(mBOLD = mean(BOLD,na.rm = T),
            seBOLD = sd(BOLD,na.rm = T)/ sqrt(n()))
## `summarise()` has grouped output by 'subject', 'voi', 'method', 'tr'. You can
## override using the `.groups` argument.
ggplot(nacc_tc_b4_summary_cueval %>% 
         filter(voi == "nacc8mm",
                cue_value %in% c("+$5","+$0")), aes(x = tr, y = mBOLD, group = cue_value, 
                                                    fill = cue_value))+
  geom_line(aes(color = cue_value))+
  geom_ribbon(aes(ymin = mBOLD - seBOLD, ymax = mBOLD + seBOLD),
              alpha = 0.2)+
  facet_wrap(~subject+method, nrow = 2)+
  theme_bw()+
  scale_fill_brewer(palette = "Set2")+
  scale_color_brewer(palette = "Set2")+
  scale_x_continuous(breaks = 1:10)

nacc_tc_b4_summary_cueval <- nacc_tc_b4 %>%
  group_by(subject,voi,method,tr,cue_value) %>%
  summarise(mBOLD = mean(BOLD,na.rm = T),
            seBOLD = sd(BOLD,na.rm = T)/ sqrt(n()))
## `summarise()` has grouped output by 'subject', 'voi', 'method', 'tr'. You can
## override using the `.groups` argument.
ggplot(nacc_tc_b4_summary_cueval %>% 
         filter(voi == "nacc8mm",
                cue_value %in% c("+$5","+$0")), aes(x = tr, y = mBOLD, group = cue_value, 
                                                    fill = cue_value))+
  geom_line(aes(color = cue_value))+
  geom_ribbon(aes(ymin = mBOLD - seBOLD, ymax = mBOLD + seBOLD),
              alpha = 0.2)+
  facet_wrap(~subject+method, nrow = 2)+
  theme_bw()+
  scale_fill_brewer(palette = "Set2")+
  scale_color_brewer(palette = "Set2")+
  scale_x_continuous(breaks = 1:10)

ggplot(nacc_tc_b4_summary_cueval %>% 
         filter(voi == "nacc8mm",
                cue_value %in% c("+$5","-$5","-$0")), aes(x = tr, y = mBOLD, group = cue_value, 
                                                    fill = cue_value))+
  geom_line(aes(color = cue_value))+
  geom_ribbon(aes(ymin = mBOLD - seBOLD, ymax = mBOLD + seBOLD),
              alpha = 0.2)+
  facet_wrap(~subject+method, nrow = 2)+
  theme_bw()+
  scale_fill_brewer(palette = "Dark2")+
  scale_color_brewer(palette = "Dark2")+
  scale_x_continuous(breaks = 1:10)

# 
# View(nacc_tc_b4_summary_cueval %>% filter(tr == 4, cue_value == "+$5",
#                                           voi == "nacc8mm"))
# 
# View(nacc_tc_b4_summary_cueval %>% filter(tr == 5, cue_value == "+$5",
#                                           voi == "nacc8mm"))

nacc_tc_b4_summary_hit5 <- nacc_tc_b4 %>% 
  filter(cue_value == "+$5") %>% 
  group_by(subject,voi,method,tr,hit) %>% 
  summarise(mBOLD = mean(BOLD,na.rm = T),
            seBOLD = sd(BOLD,na.rm = T)/ sqrt(n())) %>% 
  mutate(hit = as.factor(hit))
## `summarise()` has grouped output by 'subject', 'voi', 'method', 'tr'. You can
## override using the `.groups` argument.
ggplot(nacc_tc_b4_summary_hit5 %>% 
         filter(voi == "nacc8mm"), aes(x = tr, y = mBOLD, group = hit, 
                                       fill = hit))+
  geom_line(aes(color = hit))+
  geom_ribbon(aes(ymin = mBOLD - seBOLD, ymax = mBOLD + seBOLD),
              alpha = 0.2)+
  theme_bw()+
  scale_fill_brewer(palette = "Set1")+
  scale_color_brewer(palette = "Set1")+
  scale_x_continuous(breaks = 1:10)+
  facet_wrap(~subject+method,nrow =2)+
  labs(title = "+$5, hit vs .miss")

f1 <- nacc_tc_b4 %>% 
  filter(cue_value == "+$5") %>% 
  mutate(run = ifelse(trial <= 42, "run-01","run-02")) %>% 
  group_by(subject,voi, tr,method, run) %>% 
  summarise(mBOLD = mean(BOLD,na.rm = T),
            seBOLD = sd(BOLD,na.rm = T)/sqrt(n()))
## `summarise()` has grouped output by 'subject', 'voi', 'tr', 'method'. You can
## override using the `.groups` argument.
ggplot(f1, aes(x = tr, y = mBOLD, group = method))+
  geom_line(aes(color = method),linewidth = 1)+
  # geom_ribbon(aes(ymin = mBOLD-seBOLD,ymax = mBOLD+seBOLD, fill = method),
  #               width = 0.1, alpha = 0.1)+
  facet_wrap(~subject+run+voi)+
  scale_fill_brewer(palette = "Set1")+
  scale_color_brewer(palette = "Set1")+
  theme_bw()+
  geom_hline(yintercept = 0, linetype = "dashed")+
  scale_x_continuous(breaks = 1:10)